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Optical Communications for 
Transport Aircraft 

Robert Stengel 
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The problem 

Increasing demand for radio-frequency bands from an enlarging pool 
of users (aircraft, ground and sea vehicles, fleet operators, traffic 
control centers, commercial radio and television) 

Desirability of providing high-bandwidth, dedicated communications 
to and from every aircraft in the National Airspace System 

Need to support communications, navigation, and surveillance for a 
growing number of aircraft 

Improved meteorological observations by use of probe aircraft 

THE SOLUTION 

Optical signal transmission support very high data rates 

Optical transmission of signals between aircraft, orbiting satellites, 
and ground stations, where unobstructed line-of-sight is available 

Conventional radio transmission of signals between aircraft and 
ground stations, where optical line-of-sight is unavailable 

Radio priority given to aircraft in weather 
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Aerospace Optical Communication 

Data communications between aircraft and ground stations could be 
supported with direct and relayed signals. Aircraft at altitude typically 
would have unobstructed line of sight to an overhead spacecraft and fre- 
quently could communicate with other aircraft at similar altitude. Fiber- 
optic links on the ground complete the data path for air-ground links 
obscured by clouds through unobscured air-satellite-ground links. 



• Opportunistic Optical Transmission 

• Distributed Network containing Free-Space and Fiber- 
Optic Links 

• Radio Transmission where Optical Link is 
Unavailable 
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Typical Cloud Cover Patterns 
Across the United States 
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Every Aircraft a Weather Probe and 
Airborne Surveillance System 

Increased data bandwidth allows greatly expanded transfer of 
information about weather conditions and individual aircraft. Observational 
data from aircraft is integrated into a real-time four-dimensional weather 
map in ground-based computers. This information, in turn, becomes 
available to all aircraft in die system. 

• D own link 

Own position and velocity vectors 

Own air temperature, pressure, and humidity 

Own wind velocity vector 

Own light intensity 

Own turbulence intensity 

Signal strengths from electrical activity and beacons 
Airborne hazard status monitoring and alerts 
Desired alternate flight plans 

• Uplink 

Air temperature, pressure, and humidity fields 
Wind and turbulence fields 
Cloud cover 
Traffic alerts 

Ground/satellite-based hazard status monitoring and 
alerts 

Arbitrated alternate flight plans 
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Research Issues 


Numerous technical, operational, and institutional issues must be 
resolved before the suitability of optical communications for aircraft can be 
fully assessed. Many of these are topics for basic and applied research. 

• Optical signal generation, transmission, and detection 

• Coherence, filtering, power, multiplexing, and coding 

• Coupling between optics and electronics 

• Communication coverage modeling 

• Telescope field of view, pointing, acquisition, and 

tracking 

• Free- space/fiber-optic networking and data-relay 
protocols 

• Architectures for CNS and ATM 

• Interfaces with related systems 

• Integration within an Intelligent Aircraft/Airspace 

System 
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Robustness of Solutions to a Benchmark Control Problem 


Robert F. Stengel* and Christopher I. Marrisont 
Princeton University , Princeton, New Jersey 08544 


The robustness of 10 solutions to a benchmark control design problem presented at the 1990 American Control 
Conference has been evaluated. The 10 controllers have second- to eighth-order transfer functions and have been 
designed using several different methods, including //<* optimization, loop-transfer recovery, Irffaglnary-axls 
shifting, constrained optimization, structured covariance, game theory, and the Internal model principle. 
Stochastic robustness analysis quantifies the controllers' stability and performance robustness *lth structured 
uncertainties in up to si* system parameters. The analysis provides Insights Into system response that are not 
readily derived from other robustness criteria and provides a common ground forjudging controllers produced 
by alternative methods. One Important conclusion Is that gain and phase margins are not reliable Indicators of 
the probability of Instability. Furthermore, parameter variations actually may Improve the likelihood of achiev- 
ing selected performance metrics, as demonstrated by results for the probability of settling-lime exceedance. 


Introduction 

C ONTROL systems should be designed to maintain satis- 
factory stability and performance characteristics not only 
at nominal operating points but over a range of parameters 
that encompasses system uncertainty. These systems should be 
robust, but there is a limit. Unbounded robustness is no more 
attractive than inadequate robustness, because nominal per- 
formance and insensitivity to parameter variations tend to 
produce conflicting design requirements. Hence, the degree or 
robustness that must be furnished for satisfactory operation is 
related to the system variations that are most likely to occur. 

Measures of robustness should be easily understood and 
should be directly connected to control design objectives. They 
should be consistent with what is known about the structure 
and parameters of the plant’s dynamic model. 1 hesc goals are 
best served when robustness is expressed in terms of the like- 
lihood that commonly accepted properties fall within accept- 
able bounds and when parameter variations are expressed in 
terms of readily measured system specifications. A method of 
satisfying these evaluation criteria is presented here. 

This paper demonstrates the application of stochastic ro- 
bustness analysis (i.e., determining the probability of unsat- 
isfactory stability or performance resulting from expected 
parameter uncertainty) to solutions of the 1990 American 
Control Conference Benchmark Control Problem. 1 Stochastic 
robustness is seen to provide a useful, unifying analytical 
framework that is intuitive and has a direct, physical meaning. 

Description of the Problem 

The benchmark plant is a dual-mass/single-spring system 
with noncollocated sensor and actuator 1 ; its state-space model 
is 


y = *2 + v 

z=x 2 


( 2 ) 

(3) 


where and x 2 are the positions of the masses, xy and x 4 are 
their velocities, and u Is a control force on nt|. The plant is 
disturbed by won m 2 , and the measurement of x 2 is corrupted 
by noise v in y . The corresponding actuator and disturbance 
input/output transfer functions are 


{k/r?i\m 2 ) 

K w> = 5 2 [ 5 2 + £ { m t +Mi)/ffl|Wl] 

(I/m 2 )(5 2 + k /tn \) 

K ” 5 2 [5 2 + *('”1 + mi)/ mi m 2 


(4) 

(5) 


The plant has eigenvalues at ( + m 2 )/m , m 2 , 0,0) 
and is undamped. A single-input/single-output (SISO) con- 
troller must close its loop around 3C„ vt which has a pole-zero 
surplus of 4. The high-gain asymptote of at least one root 
lies in the right half plane for any SISO feedback compensator 
that has [ewer than two surplus zeros. Because the open-loop 
roots are on the imaginary axis, the magnitudes of root depar- 
ture angles must exceed 90 deg if marginal instability is to be 
avoided at low loop gain. 

Three design problems are posed in Ref. 1. Benchmark 
problem 1 (BP- 1) requires 1) a 1 5-s settling time for unit distur- 
bance impulse and nominal mass-spring values {m x -m 2 ^k 
- i) and 2) closed-loop stability for fixed values of mass and 
0.5 < A < 2. The second problem, BP-2, replaces the unit- 
impulse disturbance by a sinusoidal disturbance with 0.3-rad/s 
frequency but unknown amplitude and phase. Asymptotic re- 
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j eel ion of the signal should be achieved with a 20-s settling 
time Tor nominal masses and 0.5 < k <2. The third problem, 
BP-3, is like BP-1, excepl that m,, m 2 , and k are uncertain 
with mean values of I and unspecified bounds. A number of 
additional problem specifications arc left to the discretion of 
the designer. For example, it is presumed that a noise model 
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v(r) would be considered, but derails of the mode! 3 re open. 
Subjective goals include achieving reasonable performance/ 
stability robustness, minimizing controller effort, and mini- 
mizing controller complexity. 

Design Solutions and Nominal Performance 
Five papers containing design solutions appear in the Amer- 
ican Control Conference Proceedings, 2 ' 6 one paper became 
available after the conference, 2 and additional designs were 
obtained from the authors. The transfer functions for these 
controllers are presented in the Appendix. Fixed-order com- 
pensators achieving approximate loop-transfer recovery are 
motivated in Ref. 2, leading to designs A-C. An plus 
ym-axis shifting approach is taken in Ref. 3, producing design 
D. Reference 4 uses nonlinear constrained optimization to 
produce design E. Structured covariance terms are added to 
the linear quadratic Gaussian (LQG) algebraic Riccati equa- 
tions to generate design F in Ref. 5. Design G is a game-theo- 
retic controller based on linear exponential Gaussian and H 
concepts and is discussed in Ref. 6. //« controllers using the 
internal mode] principle are presented in Ref. 7 (designs H-J). 
G and J are designed to reject the sinusoidal disturbance 
(BP-2) rather than the unit impulse disturbance (BP-1). All but 
two of these designs (A and D) contain non-minimum-phase 
zeros. The benchmark criteria do not address command-input 
responses; hence, the initially reversed time response of sys- 
tems with an odd number of non-minimum-phase zeros is not 
penalized. Design G has an even number of right-half-plane 
zeros, which would not produce reversed response. 

The problem statement contains an ambiguity that could 
have affected the designers’ interpretations of satisfactory re- 
sponse. Settling lime is normally defined as an attribute of 
unit-step-function response. For example, Ogata 8 states that 
"The settling time is the time required for the response curve 
to reach and stay within a range about the final value of size 
specified by absolute percentage of the final value (usually 
2% or 5*7o)." For a second-order system the 2^o settling time 
can be precisely calculated as 4/^„, where f is the damping 
ratio and is the natural frequency of the oscillatory mode. 
However, Takahashi et al. <? found that "Exact analytical ex- 
pressions for ... settling time become prohibitively compli- 
cated for systems of order higher than two." The benchmark 
ambiguity is that the final value of a strictly stable impulse 
response (BP-1) is zero; hence, there is no steady-state value on 
which to base percentage response. 

Nominal performance characteristics of the controllers are 
summarized in Table 1, which presents compensator numera- 
tor and denominator order (Num Ord and Den Ord), two 
definitions of settling time [T$ and T**) t maximum control 
usage (w m , x ) resulting from a unit w disturbance, gain margin 
(CM), phase margin (PM), output response to 0.5/rad/s sinu- 
soidal disturbance (SR), and covariance of control response 
(^cov) to measurement noise (v) with unit standard deviation. 
All compensators are proper (the number of zeros does not 
exceed the number of poles), but three (C, D, and E) are not 


strictly proper (the number of zeros equals the number of 
poles). Hence, designs A, B, and F-J can be classified as 
low-pass filters, whereas designs C-E do not roll off at high 
frequencies. : " 

7s portrays the settling time as the time for which x 2 is 
captured within a 0.1-unit envelope about its zero steady-state 
value, given an initial unit w disturbance impulse. T ** is based 
on the damping ratio and natural frequency of the dominant 
mode and is calculated as 4/fu>„. Neither of these definitions 
adheres to the conventional definition, but each has its merits. 
T* is consistent with the BP-1 problem specification, in that it 
reflects a response to a unit w disturbance; however, it is 
amplitude dependent. T$* is independent of amplitude, but it 
is unrelated to the disturbance input and is not an accurate 
portrayal of the full system's sc tiling time in response to a unit 
step input. Table 1 indicates that only three of the compensa- 
tors satisfy a I5-s criterion by the first definition, whereas six 
c o mpeh sato r sYi avc settling times of s 2 s by the second 
definition. 

Four compensators use measurably more control than the 
others in responding to the disturbance. Increasing gain mar- 
gin generally is accompanied by increasing phase margin for 
these 10 designs (Fig. f), although the relationship is not 
monotonic. With the exception of design D, stability margins 
are less than the S-dB/30-deg rules of thumb (e.g., Ref. 10) 
often suggested as design goals for SISO systems. Sinusoidal 
disturbance rejection of most controllers is similar, although 
design D’s response is an order of magnitude smaller. Designs 
G and J, specifically intended to reject a 0.5-ra”d/s sinusoid, 
eliminate the disturbance^ompletely in the steady state. (The 
settling time in achieving this response was not evaluated.) The 
noise-response covariance oT the control is generally propor- 
tional to its peak disturbance-impulse response for strictly 
proper compensators. The three non-strictly proper compensa- 
tors have infinite control covariance in response to continuous 
white measurement noise v (with infinite bandwidth). 

Stochastic Robustness Analysis 

Stochastic robustness analysis (SRA) is based on a statistical 
portrayal of parameter variations and their effects. If parame- 
ters take a finite number of discrete values, each with known 
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Fig. 1 Nominal gain and phase relationships of the 10 controllers. 


Table 1 Nominal characteristics of 10 controllers 


Design 

Num 

Ord 

Den 

Ord 

Tj , s * 1 


u m«x 

GM, 

db 

PM, 

deg 

SR, 

db 

t/cov 

A 

2 

3 

21.0 

14.8 

0.514 

2.56 

26.7 

10.1 

6.30 

B 

2 

3 

19.5 

15.2 

0.469 

3.27 

26.8 

13.2 

13.02 

C 

2 

2 

19.7 

15.2 

0.468 

3.27 

26.5 

13.3 

OO 

D 

4 

4 

9.9 

8.8 

297.8 

15.10 

58.7 

1.47 

oo 

E 

2 

2 

18.2 

8.01 

0.884 

2.39 

22.0 

17.1 

oo 

F 

3 

4 

13.7 

22.0 

2.397 

5.15 

23.8 

13.4 

6 x 10* 

G 

5 

8 

31.3 

35.7 

1.458 

3.61 

25.4 

— oo 

173.5 

H 

3 

4 

14.9 

11.9 

0.574 

3.28 

24.5 

14.9 

2,48 

f 

3 

4 

17.8 

17.2 

0.416 

4.56 

27.5 

13.3 

0.95 

J 

5 

6 

43.2 

23.8 

1.047 

2.14 

17.5 

— OO 

77.42 


* “Defined for 0 l -unit * „ • sponse envelope for unit-impulse w. 
**“ Defined by 4/fu', {pro'.ded by B Wie). 
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or estimated probability, the analysis can be based on a finite 
number oT function evaluations, and the probabilistic result is 
exact (within the accuracy and precision of problem model* 
ing). If the parameters are continuous or the number of finite 
combinations is too large for practical computation, Monte 
Carlo evaluation can be used to estimate probabilities within 
arbitrarily small confidence intervals. If a binary judgment can 
be made ol function values (e g., satisfactory/unsatisfactory 
or stable/unstable), then the corresponding probability distri- 
bution is binomial, and confidence intervals are readily esti- 
mated from the number of function evaluations (e.g., 
Ref. II). Further details of SR A can be found in Refs. 12-17. 

Test Cases for the Benchmark Problem 

Uncertain parameters are assumed to have continuous, 
bounded, uniform, and uncorrelated probability distributions 
for this analysis. (The original problem identifies uncertain 
parameters and their bounds, making no statement about dis- 
tributions. 1 ) Three increasingly demanding sets of parameter 
uncertainties are used to test the controllers. The first two are 
specified in Ref. 1, and the third is new. 

Problem E-I: 0.5 < k <2. all oilier parameters take nomi- 

nal values, as in B P- 1 and BP-2. 

Problem E-2: 0. 5 < A: < 2 , 0.5<w t <1.5, and 0.5<m 2 

< 1.5, as in BP-3. Reference 1 does not specify limits on m , 
and values of ± 50^o are adopted here. 

Problem E-3: Same as E-2; in addition, 0<c <0. 1 , 0.9 < / 

<1.1, and 0.001 < r < 0.4 s, where c represents internal damp- 
ing between the masses; / is loop-gain uncertainty due to 
multiplicative variation in observation, control gain, or actua- 
tor effectiveness; and r is the time constant for a first-order 
lag between controller command and actuator response. Un- 
certainty in the damping ratio c increases open-loop damp- 
ing, and the time lag is always greater than the nominal value 
of zero. 

With all six parameters, the state-space model for E-3 
becomes 



x’ 

= F x 

+ G ' t/, + L 

' w 


(6) 

w here Jr ’ 

is defined as [.v, 
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L’ 

= 10 0 

0 1 /in 2 
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The compensators 

are modeled by 
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w here x t is the compensator state; u t is the actuator command; 
A , £7, C, and D are the compensator matrices; and y is .r 2 . 

Performs nef Metrics Tor the Benchmark Problem 

Robustness is best characterized by problem-dependent met- 
rics that have a direct bearing on the measurable stability and 
performance or the system. Here, they portray the likelihood 
that classical stability bbunds will be exceeded, that settling 
lime will not be achieved, and that control usage will exceed 
acceptable values. For demonstration oT SR A, parameter un- 
certainties are represented by uniform distributions within ar- 
bitrary (but reasonable) bounds. In practical application, the 


control-system designer would have similar, problem specific 
specifications to meet. 

Each of the following probabilistic performance metrics 
has a binomial distribution and is estimated using Monte Carlo 
evaluation. Uniform, bounded parameters are calculated by 
random-number generators according to the specifications of 
the previous section. The associated binomial confidence level 
depends on the number of evaluations and the value of the 
probability estimate. 13 Each estimate is the result of 20,000 
evaluations; for a probability estimate of 0.1, the 95^o confi- 
dence interval would be ± 0.004. The performance metrics arc: 

1) P t : Probability of instability. This probability portrays 

the likelihood that parameter variations force at least one 
closed-loop root into the right half plane. 

2) Pr t : Probability of settling-time exceedance. This 

probability is derived from a time-history calculation with a 
unit-impulse w input (i.e., based on T*) and estimates the 
likelihood that the actual response of z will fall outside a 
±0.1 -unit envelope after 15 s. 

3) P „ : Probability of control limit exceedance. This prob- 

ability corresponds to the requirement in Ref. 1 to minimize 
controller effort. It is the probability that peak actuator dis- 
placement will exceed a saturation limit in response to a unit 
disturbance (tv) impulse. The saturation limit was chosen to be 
one unit Tor this analysis. 

4) P } \ Probability of unsatisfactory sinusoidal distur- 
bance rejection. This probability involves the likelihood that 
the amplitude of steady-state z response exceeds one unit with 
a unit sinusoidal disturbance at 0.5 rad/s. 

Computation times indicate that current workstations are 
fast enough to execute practical SRA, and massively parallel 
computers could provide interactive turnaround. For the 
typical closed-loop system considered here, roughly 900 sets 
of eigenvalues were generated per minute per million float- 
ing-point operations per sec (MFLOP). This is drawn from 
compiled Pascal code executed on a 0.9-MFLOP Silicon 
Graphics 4D/20 workstation. The complete evaluation was 
computed at a rate of 30 sets/min/MFLOP using MATLAB 
on a Macintosh 1 1 x computer. At these rates, a 5000-MFLOP 
parallel computer (e g., 64K CM-2 Connection Machine) 
would evaluate 20,000 sets of eigenvalues in 0.25 s, and the full 
evaluation would take about three times longer. 

Results of (he Analysis 

The results of the SRA indicate a wide range of characteris- 
tics in the 10 controllers. This reflects varying emphasis in 
satisfying the problem specifications, as well as significant 
differences in compensator order and design philosophy. It 
should be emphasized that none of the controllers was de- 
signed for the express purpose of satisfying SRA criteria, and 
it is likely that each design approach could be fine-tuned to 
produce better results than those shown here. Using criteria 
that have high engineering significance, SRA provides a “level 
playing field” on which to judge the robustness of controllers 
that were designed by alternative methods. Tables 2-4 present 
results, with maximum probabilities for each evaluation prob- 
lem indicated by bold letters and minimum values underlined. 

Probability of Instability 

For the least uncertain case (E-l), over half of the con- 
trollers are estimated to have zero probability of instability, 
whereas design A has a 16^0 likelihood of instability (Table 2). 
With increasing parameter uncertainty (E-2 and E-3), all con- 
trollers have nonzero P,. The probability of design A is essen- 
tially unchanged, and design J becomes the controller most 
likely to be unstable. 

It is interesting to compare the probabilities oT instability on 
the bases of gain and phase margins, quantities often assumed 
to indicate the robustness of S1SO systems. Figures 2 and 3 
demonstrate that nominal values of GM and PM are not good 
predictors of P f . (Note that these bar charts present results for 
the 10 compensators; hence, GM and PM are not evenly dis- 


149 



Table 2 Probabilit> of instability 


Design 

E-l 

E-2 

E-3 

A 

0.160 

0.159 

0.165 

B 

0.023 

0.042 

0.039 

C 

0.021 

0.040 

0.041 

D 

0.000 

0.004 

0.059 

E 

0.000 

0.097 

0.125 

F 

0.000 

0.119 

0.224 

G 

0.000 

0.203 

0.232 

H 

0.000 

0.046 

0.099 

I 

0.000 

0.013 

0.029 

J 

0.039 

0.237 

0.245 


Table 3 Probability of 
settling-lime violation 

Design 

E-l 

E-2 

E-3 

A 

0.971 

0.962 

0.793 

B 

1.000 

0.969 

0.963 

C 

1.000 

0.968 

0.874 

D 

0.000 

0.004 

0.072 

E 

1.000 

1.000 

0.999 

F 

0.633 

0.859 

0.967 

G 

1.000 

0.999 

1.000 

H 

0.742 

0.909 

0.986 

1 

0.756 

0.918 

0.986 

J 

1.000 

1.000 

0.968 


Table 4 Probability of 
control-limit exceedance 

Design 

E-l 

E-2 

rn 

* A 

0.160 

0.159 

0.165 

B 

' £023' 

0.043 

0.047 

C 

0.021 

0.041 

0.041 

D 

T.000 

1.000 

1.000 

E 

o.ood 

0.391 

0.409 

F 

1.000 

1.000 

1.000 

G 

1.000 

0.886 

0.889 

H 

0.000 

0.133 

0.162 

1 

0.000 

0.023 

0.030 

J 

0.857 

0.542 

0.527 


tributed.) In most cases, increasing parameter uncertainty in- 
creases P f . but there are no consistent trends with GM and 
PM. Parameter variations have complex effects on the shape 
of each controller’s Nyquist plot, and these effects cannot be 
portrayed simply by changing loop gain or phase angle. 

This result brings into question the utility of transfer-func- 
tion/return-difference-matrix singular values as measures of 
the stability robustness of multi-input/multi-output (MIMO) 
systems. MIMO singular-value analysis is loosely equivalent to 
SISO gain-margin analysis (e.g., Ref. 18). Arbitrary, real pa- 
rameter variations have complicated effects on the frequency 
distributions of MIMO singular values, changing their shapes 
as well as their magnitudes. Unless the frequency distributions 
of nominal MIMO norms retain their shapes under parameter 
variation (or follow some predictable pattern), the relation- 
ships of nominal maximum or minimum values to allowable 
bounds tells little about stability robustness. Norm bounds can 
be reliably evaluated only by considering the norms of per- 
turbed systems, 

A higher compensation order does not necessarily improve 
robustness (Tables 1 and 2). The compensators with the most 
stability robustness are fourth order, and the next most robust 
controllers are second and third order. Increased nominal 
control usage, either as a consequence of a disturbance impulse 
or measurement noise, generally corresponds to decreased 
stability robustness, although design D provides a significant 
exception. 


Probability of Settling-Time Violation 

All but three of the controllers (D, F, and H) exceed the 1 5-s 
settling-time objective (defined by T$) in the nominal case 
(Table I); hence, it is not surprising that the probability of 
settling-time violation with parameter uncertainty is high as 
well (Table 3). Design D provides a notable exception; Us 
nominal T$ is 9.9 s, and P Tt is small forjtlj three evaluation 
cases. For problem E-l, half of the controllers violate the goal 
all the time, but two of the controllers with nominal T* above 
15 s (H and I) have a considerable likelihood (25foj_ of satis- 
fying the objective when the spring-constant uncertainty is 
considered. Further uncertainty (problems E-2 and E-3) re- 
duces the probability of settling-time violation for more con- 
trollers, illustrating the counterintuitive result that the effects 
of uncertainty are not always unfavorable. 

Probability of Control Limit Exceedance 

The probability of excessive control response to disturbance 
impulse P u is shown in Table 4. Over half of the nominal 
responses are within the u mtx criterion chosen for this analysis 
(Table I). Furthermore, there is an identifiable trend in the 
relationship between w ma , and P u (Fig. 4). Several controllers 
(E, H, and 1) have zero probability of violating this criterion 
for problem E-I, and designs B, C, and I retain low values of 
P u for all three problems. Designs D and F have !00^o P v in all 
three cases, which is traceable to very high nominal control 
usage. Once again, nominally marginal cases (G and J, the two 
controllers designed for rejection of the sinusoid) exhibit re- 
duced probability of exceedance for problems E-2 and E-3. 
From Eq. (n. mcYeaTed m, and m 2 decrease the effects of u 
and h\ and added damping (c) and first-order lag (t) reduce 
control peaks in some cases, reducing the probability of high 
control levels. 

■ P 51 
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Fig. 2 Probability of Instability vs gain margin for three evaluation 
problems. 
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Fig. 3 Probability of instability vs phase margin for three evaluation 
problems. 
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Fig. 4 Probability of control-limit exceedance vs nominal maximum 
control response to a disturbance impulse. 
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c) 

Fig. 5 Stochastic root locus of design II (scatter plot): a) problem 
E-l; b) problem E-2; c) problem E-3. 

Sinusoidal Response Characteristics 

When 0 dB is chosen as an upper response limit, the two 
controllers designed to reject the sinusoid (G and J) do so 
perfectly ( P f = 0), whereas all the others exceed the limit all the 
time. The transfer functions (Appendix) show that designs G 
and J effectively "notch” the 0.5-rad/s disturbance-input fre- 
quency to produce these results. Without notch filters the re- 


maining controllers cannot give special attention to discrete- 
frequency inputs, and their frequency response of -0.5 rad/s 
always exceeds 0 dB. IT the frequency of the sinusoidal dis- 
turbance were uncertain, the notch fillers could be less effec- 
tive, but there would be little change in the response of the 
other controllers. 

Stochastic Root Loci and Parametric Histograms 

Graphical results give insight into the nature and causes of 
possible instability. The stochastic root locus is an 5 -plane plot 
of the eigenvalues that result Trom each Monte Carlo evalua- 
tion, expressed either as a two-dimensional scatter plot of 
closed-loop roots or an oblique three-dimensional view of the 
density of roots within subspaces oF the s plane. 13 The former 
plot is easily generated from the calculations, and the latter has 
the advantage of showing the distribution along the real axis. 17 
In addition, histograms of the parameters associated with in- 
stability can be related to origins of the problem. 

Scatter plots for design H show the progression of eigen- 
value uncertainty from problem E-l to E-3 (Fig. 5). For prob- 



lig. 6 Stochastic root locus of design II (three-dimensional view): 
a) problem E-l; b> problem E-2: c) problem E-3. 
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lem E-l only a single parameter varies (the spring constant k). 
The distribution follows the conventional root locus (with 
nominal closed-loop locations indicated by x), although the 
density of roots varies along the curves. The pairs of roots near 
the origin are most closely associated with the plant, whereas 
the higher-frequency roots are compensator modes. None of 
the root loci extend into the right half plane, and Pj is zero 
(Table 2). Three parameters vary in problem E-2, and the 
stochastic root locus becomes an areal distribution of roots, 
some of which extend into the right half plane (Fig. 5b). 
Because the parameter variations are bounded, there are crisp 
edges to the distributions. The unstable cusps at 0.6 and 2.6 
rad/s can be associated with plant and controller modes. Fur- 
ther parametric uncertainty (problem E-3) broadens the distri- 
butions and increases the probability of instability. 

The same information is presented in unsmoothed three- 
dimensional Form in Fig. 6 (upper half plane only), which 
shows the distribution of real roots as well. The three-dimen- 
sional representation is especially effective when displayed on 
a graphics workstation that allows the viewpoint to "fly 
around” the distribution. 

To see which parameter values are associated with instabil- 
ity, the values are recorded whenever the system is found to be 
unstable. These values are collected in intervals, the number of 
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Fig. 7 Parameter histograms For all unstable cases, design H: 

a) problem E-l; b) problem E-3, 



Fig. 8 Parameter histograms for high-frequency unstable cases, de- 
sign H, problem E-3. 



values in each interval is counted, and the resulting histogram 
provides an estimate of the conditional probability density 
function for each parameter, ff a parameter has little effect on 
stability, then the histogram should show the same distribution 
as produced by the random number generator — in this case, a 
uniform distribution. If particular values of the parameter 
increase the probability of instability, the histogram has higher 
values in that region. 

For design H and problem E-2, instability often occurs when 
the masses have low values but never occurs with high values 
(Fig. 7a). Low mass values increased the probability of insta- 
bility for all the designs. Extreme values of the spring constant 
also are associated with instability, low values having the edge 
in this example. 

For problem E-3 (Fig. 7b), the distributions become less 
crisp, as otherwise unstable values of mass can be stabilized by 
damping and otherwise stable values of mass can be desta- 
bilized by increased loop gain or first-order lag. The spring 
constant shows a slight bimodal distribution due to the two 
modes of instability with roots of approximately 0.6 or 2.6 
rad/s. This can be seen by recording the parameter values only 
when the system is found to be unstable and the unstable roots 
have a high frequency. The resulting histograms (Fig. 8) show 
that there are unstable high-frequency roots only if the spring 
constant is high and the damping is low. With increased damp- 
ing, there is no high-frequency instability. 

These results can be used in three ways. The probability of 
instability could be reduced if it were possible to ensure that 
the plant parameters did not move into the areas that are 
found to cause problems. This might be the result of improved 
quality assurance on the important parameters or by shifting 
the mean of the parameter variation. If it is not possible to 
affect the actual parameter variations, then the control system 
could be redesigned using the problematic values of parame- 
ters as nominal values. For example, the control system could 
be redesigned using nominal values of 0.7 for the masses. 

A third use of the distributions c§n occur if one of the 
varying parameters represents a control design parameter. For 
instance, if the loop gain / were treated as a design vari- 
able, then it is clear that attenuating the gain would reduce the 
probability of instability. This alternative is demonstrated 
using design D. It has been seen that design D had generally 
good robustness but very high actuator use. Peak actuator 
usage can be reduced by reducing the loop gain, and the effect 
of gain attenuation on robustness subject to problem E-2 is 
shown in Fig. 9. For this analysis, only 100 Monte Carlo 
evaluations were carried out per design point, but the results 
show clear trends. As the gain is reduced, the probability of 
control saturation is reduced without significant increase in 
Pi or P Ti until the attenuation reaches 0.6, when P Ti begins 
to increase. Reducing the gain further produces a clear trade- 
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off between the probabilities of control saturation and set- 
tling-time violation. References 19 and 20 present similar 
methods of control system design based on search and statisti- 
cal evaluation. 

Conclusions 

Stochastic robustness analysis of 10 controllers designed for 
the ACC Benchmark Control Problem provides useful quan- 
tification of stability and performance sensitivities to parame- 
ter variations. The SRA method is flexible and can be tailored 
to the design requirements and system specifications of partic- 
ular control problems. Qualitative selection of the best con- 
troller depends on the relative importance oT several metrics, 
which are readily described in a probabilistic framework. 

Several conclusions can be drawn from this analysis. The 
analysis shows that gain and phase margins are not good pre- 
dictors of the relative stability robustness of different 51SO 
controllers, because robustness is tied closely to the actual 


plant uncertainties and their effects on (implied) Nyquist con- 
tours. This result implies that robustness analyses based on 
singular-value analysis of MIMO systems may have similar 
limitations. Nominal settling time did not give a good indica- 
tion of the likelihood of exceeding settling time limit, princi- 
pally because most nominal values already exceeded the limit. 
Although this result may be an artifact of the settling-time 
definition (T/), it reveals the counterintuitive result that 
uncertainty may improve the probability oT remaining within 
a predefined limit. The relationship between maximum control 
response to a disturbance impulse and the probability of ex- 
ceeding a control limit is more direct, as most nominal values 
were about half the limit value. Stochastic root loci and pa- 
rameter histograms provide insight about the likely positions 
of the closed-loop roots and the parameter variations that lead 
to instability, and they suggest ways of improving plant and 
controller design. 


Appendix: Transfer Functions of the Ten Compensators 


Design A: 


Design B: 


Design C: 


Design D: 


Design E: 


Design F 


40.42(5 + 2.388)(s + 0.350) 

(5 + I63.77)[s j + 2(0. 50l)(0. 924)5 + (0.924) 2 J 


42.78(5 - l.306)(5 +0.1988) 

(5 + 73. 073) [ 5 2 + 2(0.502)(I.I82)5 +(1.182) I J 


0.599(5 - 1 .253)(5 +0.1988) 

[s 2 + 2(0.502)(l . 182)5 + (1 .182)-] 

19881(5 + I00)(5 +0.2I2)[5 2 + 2(0.I73)(0.733)5 4 (0.733) 2 ] 

[5 2 + 2(0.997)(5!.16)5 + (5 1.16) ! j[5 2 + 2(0.8381(16.44)5 + (I6.44)-’| 

S.369(5 -0.348)(5 +0.0929) 

[5 2 + 2(0.832X2.21)5 + (2.2 l) 2 J 

2246.3(5 + 0.237)[5 2 - 2(0.32)(l .064)5 + (1 .064) 2 ] 

(5 + 33. 19)(5 + 1 1 .79) [5 2 + 2(0.90X2.75)5 + (2.75) 2 ] 


Design G: 

4430(5 + 0.08)(5 - 0.44)(5 - 2.83)[5 2 - 2(0. 102)(0.49)5 + (0.49) 2 ] 

[[ 5 2 + 2(0.70)(H.17)5 + (II. 17) 2 ][5 2 + 2(0.89)(3. 67)5 + (3. 67) 2 j(5 2 + 2(0.29X3. 1 1)5 + (3.1 l) 2 ][5 2 + (0.5) 2 ]j 

Design H: 

2.13(5 +0.145X5 - 0.98)(5 + 3.43) 

[5 2 + 2(0.82)(1 .59)5 + (1 ,59) ! ] [ 5 ' + 2(0.46)(2.24)5 + (2.24)’] 

Design 1: 

16.1(5 +0.134X5- I.I74)(5 + 1.46) 

[5 ! + 2(0.82X1.05)5 + ( 1 .05) 2 J [j 1 + 2(0. 5)(2. 18)5 + (2. 18) 2 j 


Design J: 

51.47(5 +0.06)(5 -0.21 )(5 + 5.4i) [5 2 - 2(0.07)(0.5 1)5 + (0.5I) 2 } 

|5 2 + 2(0.72X2.05)5 + (2.05) J ] [ 5 2 + 2(0.68X5.21)5 + (5.21) 2 ][5 2 + (0.5) 2 j 
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Stochastic Prediction Techniques for 
Wind Shear Hazard Assessment 




D. Alexander Stratton* and Robert F. Stengelf 
Princeton University, Princeton, New Jersey 08540 


The threat of low-altitude wind shear has prompted development of aircraft-based sensors that measure winds 
directly on an aircraft’s intended flight path. Measurements from these devices are subject to turbulence inputs 
and measurement error, as well as to the underlying wind profile. In this paper stochastic estimators are 
developed to process onboard Doppler sensor measurements, producing optimal estimates of the winds. A 
stochastic prediction technique determines the level of aircraft energy performance from the wind estimates. 
Aircraft performance degradation algorithms presented are based on optimal estimation techniques. The predic- 
tion algorithm must balance wind shear detection performance and turbulence rejection capability, as illustrated 
in simulations of microbursl wind shear and severe turbulence environments. 


Introduction 

S TRONG variable winds in the airport vicinity can cause 
unacceptable deviation of aircraft from their intended 
flight path. Known as low-altitude wind shear, this threat has 
caused at least 24 aviation accidents in the last 25 years. 1 
Efforts to promote the avoidance of severe wind shear have 
focused on improving flight crew training programs, 2 under- 
standing the meteorology of wind shear, 3 ' 5 and developing 
technology to detect wind shear in the terminal area. Ground- 
based sensor systems to measure airport-vicinity winds are 
being developed and installed at major airports, 6 7 along with 
techniques to automatically identify a wind shear and predict 
its formation. 810 Sensors to detect wind-shear-induced flight- 
path deviations are being installed on aircraft, 11 12 and for- 
ward-looking sensors to detect wind shear in front of the air- 
craft also are under development. 1315 Interpretation of this 
information in the cockpit is a topic of current research. 

As the amount of available information grows, accurate 
interpretation of the information by flight crews becomes 
more challenging, particularly during periods of high work- 
load. Artificial intelligence technology provides a basis for a 
cockpit aid to assist flight crews in avoiding low-altitude wind 
shear. An expert system, the Wind Shear Safety Advisor, 16 
depicted schematically in Fig. 1, will operate in real time, 
accepting evidence from onboard and ground-based sources, 
perhaps facilitated by a direct data link (represented by a dot- 
ted line in Fig. 1). The goal of this system is to increase flight 
crew situation awareness and decision reliability by summariz- 
ing information from a variety of information sources. 

In the absence of direct measurements of the winds, a deci- 
sion to avoid wind shear must be based on discrete alerts from 
wind shear detection systems and meteorological evidence. 
Various levels of reliability associated with this indirect evi- 
dence complicate the risk assessment process. A probabilistic 
model of this process has been developed that incorporates 
statistics from meteorological studies and reliability statistics 
for wind-shear-alerting systems. 17 The model can manage the 
uncertainty associated with indirect evidence, providing mean- 
ingful estimates of risk. 
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If onboard measurements of the winds were available, a 
hazardous level of wind shear could be identified by deter- 
mining whether the level of some hazard metric, based on the 
wind measurements, exceeds a threshold. Hazard metrics 
considered previously include maximum horizontal winds 3 
and F-factor, 14 which relates wind shear to aircraft perfor- 
mance. Computation of the hazard level is complicated by 
uncertainty surrounding the wind measurements, including 
turbulence and measurement errors. In this paper Kalman fil- 
ters are developed to produce optimal wind estimates from 
onboard wind sensors, based on a stochastic wind model. 
These algorithms are demonstrated in a simulated microburst 
wind shear environment. 

From the wind estimates, predictions of the aircraft’s per- 
formance degradation can be made using stochastic predic- 
tion techniques. 18,19 In addition to the predictions themselves, 
these techniques produce measures of the possible error in the 
predictions due to turbulence and limitations of the measure- 
ment devices. In this paper a Kalman-filter-based prediction 
technique to predict F-factor and aircraft performance degra- 
dation is demonstrated in simulated microburst wind shear 
encounter. The response characteristics of the prediction tech- 
nique must provide significant response to severe wind shear 
and limited response to turbulence. In this paper stochastic 
prediction techniques with different design parameters are 
demonstrated in a simulated microburst wind shear and severe 
turbulence environments. 

Probabilislic Reasoning in Artificial Intelligence 

The power of an intelligent system rests in its ability to 
produce meaningful conclusions by reasoning, i.e., by apply- 
ing knowledge stored in the system to available evidence, in 
probabilistic models of reasoning, knowledge is stored in the 
form of probabilities, and Bayes’s rule 20 and the axioms of 
probability 21 are used to condition these probabilities on evi- 
dence. When several pieces of evidence are supplied, the appli- 
cation of Bayes’s rule is complicated by dependencies between 
pieces of evidence. A structure to these dependencies must be 
provided for efficient reasoning. In Bayesian network repre- 
sentation 22 a graphical representation provides this structure, 
such as the one for wind shear avoidance graphed in Fig. 2. 
Nodes in the diagram represent discrete random variables, and 
the links between them represent sets of conditional probabil- 
ities used during reasoning. The network representation ena- 
bles efficient probabilistic reasoning because all of the depen- 
dencies between variables are specified by the links. 

The network of Fig. 2 was developed using guidelines for 
wind shear avoidance presented in the FAA’s Windshear 
Training Aid document, 2 which was written by a team from 
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Fig. 1 Wind shear safe*) advisor schematic diagram. 



Fig. 2 Graphical representation of a Bavesian network for wind 
shear avoidance. 


the airframe industry with the support of airlines, the govern- 
ment, and academia. The network model incorporates statisti- 
cal results from the NIMROD, 3 JAWS,- 4 and FLOWS 5 stud- 
ies and for the enhanced Low-Level Windshear Alert System 
(LLWAS) evaluation. 7 Demonstrations of the network 17 show 
that it can approximate the subjective judgments required to 
establish the possible presence of wind shear. 

A probabilistic model establishes a scientific basis for 
the Windshear Training Aid avoidance guidelines. Since the 
completion of the Windshear Training Aid, a variety of new 
ground-based and airborne wind shear detection systems are 
being devleoped, such as the Terminal Doppler Weather Radar 
(TDWR) system. The probabilistic model can be expanded to 
include statistics from new detection systems established dur- 
ing their ev aluation. New knowledge gained from meteorolog- 
ical studies, such as geographical variation of wind shear fre- 
quency, can also be included. 

Kalman Filler Development 
for Doppler Wind Measurements 

Airborne sensor technology with the capability to detect 
wind shear in front of the aircraft is currently under devel- 
opment, including Doppler radar, 13 Doppler lidar, 14 and in- 
frared 15 technology. Doppler devices measure a shift in fre- 
quency of radar or light waves emitted along a radial line, 
measuring the component of wind velocity parallel to that line. 
Operational devices could provide measurements of head 
winds or tail winds at a series of locations along the aircraft’s 
intended approach or takeoff path. For example, airborne 
Doppler radars could provide measurements spaced at - 500- 
ft intervals over a range of 3-5 miles, spanning 50-100 s of 
flight at approach speed. 13 This sequence of measurements 
contains the effect of turbulence and is corrupted by measure- 
ment noise as well. A bank of Kalman filters can improve the 
accuracy of hazard estimates based on successive measurement 
sequences, minimizing measurement noise and accounting for 
correlation in the wind field using a stochastic model. 

As the aircraft travels down the flight pat- measurements in 
successive sequences are offset by a distance d (F ig. 3), which 


is assumed to be small relative to the distance between adjacent 
range gates L. At a given time, a sequence of measurements is 
obtained. Each member of this sequence represents the aver- 
age value of the radial wind component in an interval of length 
L at that time. 

A first-order Markov model for the turbulent winds can be 
based on the Dryden power spectrum for horizontal turbu- 
lence, given by Ref. 23 as 

. , , (2Ly\ i 

JlTTiZ^I «» 

Parameters of this model include the turbulence scale length L u 
and the root-mean-square turbulence amplitude a u . The corre* 
spending discrete Markov sequence is 

“'r, = exp( - d u + V I - exp( - 2d j jj* . , (2) 

where d u is (he ratio of d to L „ . The ij is a normally-distributed 
white noise sequence with mean and variance: 

£■(»?*)= 0 (3) 

EWk\=ol/* (4) 

This mode! uses the discrete white noise sequence tj to approx- 
imate the integrated effect of continuous white noise. Figure 4 
presents the autocovariance function associated with Eq. (1), 
along with the autocovariance function of the sequence of 
Eq. (2), indicating the agreement of the turbulence models. 

With the assumption that measurement noise is super- 
imposed on the radial wind components, the measurement at 
range gatey during measurement sequence k , z Jk can be related 
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Fig. 4 Comparison of Dryden turbulence spectrum autocovartance 
function and autocovariance function of discrete turbulence model. 
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to the corresponding scalar radial wind component by the 
relationship 

Zjk = uv y4 + n jk + v, k (5) 

This can be rewritten as 

Zjk - W r jk + fijk ( 6 ) 

where V f is the aircraft’s inertial speed at the time of measure- 
ment sequence k , and i }k has this bias subtracted out. Error in 
the inertial speed estimate, n i( , which is made from onboard 
measurements, is added to n jk to produce : 

(7) 


takes the same form as Eqs. (8-12), except that the distance 
between range gates L is used as the distance between measure- 
ments d . 


Hazard Metrics and Stochastic Prediction 

The detection of the presence of a wind shear can be based 
on the output of the stochastic estimators. A reasonable ap- 
proach to detecting wind shear is to predict whether the level 
of some hazard metric based on the w ind estimates will exceed 
a threshold. The F-factor hazard metric relates wind shear to 
aircraft air-referenced specific energy rate, which is defined by 



The measurement error n is assumed to be a zero-mean, nor- 
mally distributed white noise sequence, with a known constant 
standard deviation <r„. 

With the aforementioned assumptions, an estimator dedi- 
cated to each range gale can be constructed in the form of a 
Kalman filter. From the measurement i jk , each Kalman filter 
constructs an estimate *V, 4 ( + ) and a variance /?,*( + ), which 
is a measure of the uncertainty in iv, ( + ), in three steps. First, 
the state estimate and variance form the previous measure- 
ment sequence, . ,( + ) and p ]k _ d + ), are extrapolated ac- 
cording to 


<%(-) = exp(-rf u )w v ( ( + ) (8) 


where V 0 is the airspeed, h is aircraft altitude, and g is the 
gravitational constant. Using longitudinal aircraft equations 
of motion and assuming small flight-path angles, it can be 
shown 14 that 


d£, (T-D)V a 

-<o-— -*(')•'. 


(H) 


where T is thrust, D is drag, and W is aircraft weight. ${t) is 
the F-factor, defined as 
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(15) 


Pjd - ) « exp( - 2d u )p jk .,( + )+ [l - exp( - 2d u )\ o* /ir (9) 

Equation (8) is obtained by taking the expected value of Eq. 
(2). Note that Eq. (9) is an approximation to the integrated 
effects of continuous white noise. Next, the extrapolated vari- 
ance Pjk{~) is used to compute a gain K jk : 


K jk 


Pjk(-) 

Pjk{-) + o 2 „ 


( 10 ) 


Finally, the post-update wind estimate and variance are com- 
puted: 


*% ( + )-<%(-) + K jt ft* - ( - )] (11) 

P,*i + )= [p,*(-) + «2J/1/»,.(-)<»21 (12) 

The Kalman filters compute a weighted average of the wind 
measurements obtained at each range gate, compensating for 
the movement of the sensor platform by making an assump- 
tion of frozen Dryden turbulence in the interval between the 
measurements. Wind shear estimates are updated at each mea- 
surement step, compensating for turbulence and weighing cur- 
rent and prior information according to its relative uncer- 
tainty. Because each range gate’s state estimator h decoupled 
from the others, the computation could be performed on a set 
of identical processors running in parallel. This decoupling is 
achieved as a consequence of the Markov property of the wind 
model: the probability distribution at a given u'ind state w, 
is conditionally independent of w, J# given the closer state 
w r jk _ | • This assumption could be relaxed, coupling adjacent 
states or larger groups of states together with a corresponding 
increase in computational complexity. 

Prior state estimates and variances are required to initialize 
each filter. This may be accomplished by applying a separate 
initialization Kalman filter to the first sequence of wind 
measurements. This filter is initialized with an onboard wind 
estimate and variance at the aircraft’s location, perhaps from 
a Kalman filter processing onboard sensor measurements. 
An initial sequence of wind measurements from the forward- 
looking sensors is then processed to initialize the state and 
variance of each Kalman filter. The initialization Kalman filter 


where u\(/) is the wind component in the inertial horizontal 
direction, and w A (r) is the vertical wind component. For small 
flight-path angles, the radial wind components are approx- 
imately the same as the longitudinal horizontal wind compo- 
nents. Wind shear effects enter Eq. (14) in three ways: l) by 
changing the airspeed, 2) by altering the drag, and 3) directly 
through For conditions typical of jet transport flight 

through severe wind shear, only the direct impact of T(r) is 
significant. Prediction of aircraft specific energy along the 
intended trajectory appears to involve the prediction of air- 
speed, but using a constant nominal value of airspeed in Eq. 
(15) introduces a small, conservative error. 

The first component of $ in Eq. (15) is proportional to the 
rate of change of the horizontal wind component. If the wind 
field is assumed stationary, prediction of 5 along the intended 
trajectory could be made by differencing adjacent wind esti- 
mates: 


5 , = l/£(*v,, (16) 

This would amplify high-frequency noise, resulting in exces- 
sive prediction error. Alternatively, predicted energy deviation 
and J can be computed by a Kalman filter algorithm using the 
wind estimates as inputs. J is obtained through a weighted sum 
of the radial wind estimates, with the weights selected by defi- 
nition and minimization of a suitable cost function. 

An important limitation of Doppler wind measurement de- 
vices is their inability to measure winds perpendicular to the 
direction of the Doppler pulse. As a consequence, the second 
component of T in Eq. (15), due to vertical winds, is not 
measured by the device. In downburst wind shears, head-tail 
wind shear is produced by vertically descending winds that 
flow outward as they near the ground. These downdraft winds 
pose a hazard to the aircraft that the Doppler sensors cannot 
directly measure. Current research is attempting to model the 
vertical wind as a function of the horizontal wind for hazard 
estimation. In the simple downburst model of Ref. 23, the 
correlation between horizontal and vertical winds depends on 
the size of the downdraft, the altitude, and the distance from 
the downburst core. In a well-measured and well-studied mi- 
croburst, four major downdraft regions were found. 24 As the 
relationship between horizontal and vertical winds remains to 
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be established, the present study is based on radial wind alone. 
If a consistent correlation between vertical wind and radial- 
wind measurement is found, vertical wind could be added to 
the stochastic model. 

To predict the wind-shear-induced energy deviation £ ih , 
Eq. (14) can be integrated across a typical range gate j % result- 
ing in the recursive form 



and 

EK)=< (25) 

With the previously given model, prediction of the hazard 
level can be made from the output of the estimation Kalman 
filters after each measurement sequence. The wind estimates 
are processed using a recursive procedure based on the Kalman 
filler. 18 19 The prediction is initialized with onboard estimates 
of w x and T 0 - Predictions of £* JH and denoted £ su and 
are made for each range gate using the recursive equations 


where V, is average inertial speed of the aircraft. is modeled 
as a stationary process driven by a discrete random sequence: 

= + 1,-1 0 g ) 

where rj is a normally distributed white noise sequence with 
zero mean and standard deviation a n . This standard deviation 
is a design parameter that alters the response characteristics of 
the prediction filter, as demonstrated by simulation. Equa- 
tions (17) and (18) may be written in vector-matrix form: 



where 

< 20 > 

The relationship between prediction and estimation is obtained 
by substitution of Eq. (15) into Eq. (14) and integration from 
the aircraft (denoted with subscript 0) to a typical range gate 
j . This results in the equation 

~ 0 = ^ (£**, - ^ ** J + ( "f) H * *' ^ ^ 

If the prediction is initialized with the condition 

£,h 0 = (22) 


then Eq. (21) may be rewritten as 



(23) 


In this paper vertical wind is modeled as a normally dis- 
tributed white random sequence, uncorrelated with the radial 
winds, with mean and variance 


*KJ-° 


(24) 
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(27) 


These equations involve two gains, K Ej and K that are com- 
puted at each step based on the covariance propagation and 
filter gain computations of the Kalman filter. ,,J * The design 
parameter o n influences the size of these gains, influencing the 
response characteristics of the prediction filters. 


Simulation of Stochastic Prediction Techniques 
The stochastic estimation and prediction algorithms are 
demonstrated using a batch simulation of aircraft encounters 
with downbursi wind shear and with severe turbulence. For 
each simulation, two different predictions are made, based on 
different choices of the design parameter cr,. The wind shear is 
modeled by the Oseguera-Bowles stagnation-point-flow down- 
burst model, 25 and severe turbulence is modeled using the Dry- 
den spectrum as presented in Ref. 26. A twin-jet transport 
aircraft is represented by a point-mass longitudinal model, 27 
trimmed along an approach path at a constant airspeed of 160 
Kts. Normally distributed white noise is superimposed on mea- 
surements to simulate Doppler sensor error. Table 1 lists the 
parameters of the simulation. 

The wind shear simulation is initiated with the microburst 
just out of the sensor's detection range. Figure 5 depicts the 



Table 1 Simulation parameters 


Aircraft initial conditions 

Airspeed, V a 

160 Kt 

Altitude, h 

2000 ft 

Inertial flight-path angle. 

-3 deg 

Distance io microbursi core 

20.100 ft 

Doppler sensor 

Range gate separation, i 

500 ft 

Distance between sequences, d 

27 ft 

Noise standard deviation, a n 

1 ft/s 

Distance to aircraft 

20,000 ft 

Turbulence 

rms turbulence intensity, o » 

2.7 ft/s 

Turbulence scale length, L u 

1000 ft 

Microburst 

Downdraft radius 

2070 ft 

Maximum horizontal winds 

- S.4 ft/s 

Height of boundary laser 

131 ft 


Fig. 5 Comparison of microburst model headwind-tailwind compo- 
nent of F-factor with predicted F-faclor. 



Fig. 6 Comparison of aircraft energy deviation due to headwind-tail- 
wind shear and predicted energ> deviation. 
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Fig. 7 Comparison of F-Faclor prediclions in severe Dryden turbu- 
lence. 


situation 10 s later, comparing the predicted hazard metric T av 
along the flight path with the model’s £F a , component due to 
the headwind/tailwind shear alone. The predictions agree well 
with the model’s head/tail wind component of T avt but the 
peak magnitude of the prediction is attenuated due to the finite 
bandwidth of the prediction algorithm. In addition, the dis- 
tance between the aircraft and the wind shear is overpredicted 
due to phase shifting. With a lower value of o,, the estimators 
have lower gains, and these effects are more pronounced. If 
a wind shear warning were issued each time a critical value of 
was exceeded, the algorithm with higher o n would have a 
greater chance of positively identifying severe wind shear. 

For the same simulation, Fig. 6 compares the predicted en- 
ergy deviation, normalized as an airspeed deviation, and the 
energy deviation due to the component of the wind shear. 
Although the error in prediction of distance to the microburst 
is greater for the lower value of a,, both predictions perform 
favorably in predicting peak energy loss. However, the total 
energy loss to the aircraft is greater than either prediction, due 
to the effect of the unobserved downdraft winds. 

Figure 7 compares the predicted hazard metric 3 av for each 
of the prediction designs in severe Dryden turbulence. The 
higher choice of a, results in greater response to turbulence. 
If wind shear warnings were issued each time a critical value 
of fF av was predicted, the algorithm with higher a, would issue 
more frequent false alarms. The optimization of a prediction 
algorithm must take into account both detection performance 
and false alarm prevention. Wavelengths corresponding to 
severe wind shear should be passed, but short wavelength 
disturbances that do not affect the flight path should be 
eliminated. 

Conclusions 

Doppler wind sensors can provide advance warning of a 
wind shear threat, but wind measurements are influenced by 
turbulence and measurement error. Optimal estimation pro- 
vides a framework for minimizing the error of wind estimates 
given a hypothesis of the wind field structure. The estimation 
procedures presented here assume a structure to the local wind 
field at each range gate of the Doppler sensor, resulting in a 
bank of parallel Kalman filters. A-first-order Markov turbu- 
lence model accounts for spatial correlation in the wind field 
due to turbulence. Measures of uncertainty are produced dur- 
ing the optimal estimation process. Stochastic prediction tech- 
niques are used to predict (he impact or estimated winds on the 
energy performance of the aircraft. These techniques extend 
naturally to multiple Doppler sensors and could be expanded 
to predict other quantities such as altitude deviation error and 
touchdown dispersion error, given a nominal model of pilot 
compensation. 

If wind shear warning is based on a critical threshold value 
of a hazard prediction, the detection reliability depends on 
the design of the prediction algorithm. Kalman-fiiter-based 
designs may be band limited, identifying areas with a sus- 
tained level of substantial wind shear. To further refine the 


algorithm, a ^mparative analysis of prediction algorithm de- 
signs can be conducted, using an ensemble of representative 
severe wind shear models. The potential for false warning in 
severe turbulence also can be compared. Both threshold and 
design bandwidth may be chosen to further optimize detection 
reliability. 

Hazard prediction from Doppler sensors can provide the 
sole basis for a wind shear alert, but the lack of vertical wind 
estimates limits the alert’s reliability. Other sources of infor- 
mation could improve the reliability of Doppler-based stochas- 
tic predictions through adaptive prediction techniques. More- 
over, threshold exceedance of a hazard prediction could be 
viewed as uncertain evidence supporting a hypothesis of severe 
wind shear in the Bayesian network. With the reliability of 
threshold exceedance as evidence established through statisti- 
cal analysis, hazard prediction can be incorporated into a 
probability-based expert system for wind shear avoidance. 
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Stochastic robustness, a simple technique used to estimate the robustness or linear, time-invariant systems, is 
applied to a twin-jet transport aircraft control system. Concepts behind stochastic stability robostness are 
extended to stochastic performance robustness. Stochastic performance robustness measures based on classical 
design specifications and measures specific to aircraft handling qualities are introduced. Confidence intervals for 
comparing two control system designs are presented. The application of stochastic performance robustness, the 
nse of confidence intervals, and tradeoffs between performance objectives are demonstrated by means of the 
twin-jet aircraft example. 


Introduction 

S TANDARD linear control system design techniques rely 
on accurate models of the system to be controlled. Be- 
cause models are never perfect, robustness analysis is neces- 
sary to determine the possibility of instability or inadequate 
performance in the face of uncertainty. Robustness to these 
uncertainties, parametric or unstructured, is normally treated 
deterministically and often without regard for possible physi- 
cal variations in the system. Consequently, ovcrconservative 
control system designs or designs that are insufficiently robust 
in the face of real-world uncertainties are a danger. 

Stochastic robustness analysis (SRA), a simple technique to 
determine the robustness of linear, time-invariant systems by 
Monte Carlo methods, was introduced in Ref. 1 and presented 
in detail in Refs. 2 and 3. These references described stochastic 
stability robustness analysis and introduced the probability of 
instability as a scalar measure of stability robustness. Confi- 
dence intervals for the scalar probability of instability were 
presented, and the stochastic root locus, or probability density 
of the closed-loop eigenvalues, graphically portrayed robust- 
ness properties. Because it uses knowledge of the statistics of 
parameter variations directly, SRA provides an inherently pre- 
cise yet simple characterization of robustness. The physical 
meaning behind the probability of instability is apparent, and 
overconservative or insufficiently robust designs can be avoided. 
Applications of SRA to full-state feedback aircraft control 
systems were described in Ref. 4. The results presented there 
illustrated the use of stochastic stability robustness techniques 
in comparing control system designs and in including finite-di- 
mensional uncertain dynamics. 

Concepts behind stochastic stability robustness can be ex- 
tended to provide insight about control system design for 
performance. Design specifications such as rise time, over- 
shoot, settling time, dead time, and steady-state error nor- 
mally are used as indicators of adequate performance and lend 
themselves to the same kind of analysis as already described. 
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Concepts of stochastic stability robustness analysis can be 
applied to these criteria giving probabilistic bounds on scalar 
performance criteria. Metrics resulting from SRA can be re- 
lated to controller design parameters, thus providing a foun- 
dation for design tradeoffs and optimization. Extensions and 
uses of stochastic performance robustness in aircraft control 
system design and analysis are described in the following, and 
they are illustrated by means of an example. 

Stochastic Performance Robustness 

Stochastic stability robustness analysis is based on Monte 
Carlo analysis of the probability of instability P , and associ- 
ated confidence intervals, given a statistical description of pa- 
rameter uncertainty. 2 ' 4 Because the stability test is binomial 
(i.e., the outcome of each Monte Carlo evaluation takes one 
of two values: stable or unstable), lower L and upper U 
confidence bounds are calculated using the binomial test. 5 
While stability is an important element of robustness, perfor- 
mance robustness analysis is vital to determining whether im- 
portant design specifications are met. Adequate performance, 
such as initial condition response, command response, control 
authority, and rejection of disturbances, is difficult to de- 
scribe by a single scalar metric. Nevertheless, elements of 
stochastic stability robustness analysis apply for binomial per- 
formance metrics. 

Numerous criteria stemming from classical control concepts 
exist as measures of adequate performance. Appealing to 
these, one can begin a smooth transition from stability robust- 
ness analysis to performance robustness analysis simply by 
analyzing the degree of stability or instability rather than strict 
stability. As described in Ref. 2, one method of doing this is to 
shift the vertical discriminant line from zero to I < (or >) 0. 
Histograms and cumulative distributions for varying degrees 
of stability are readily given by the Monte Carlo estimate of 
the probability of any eigenvalue real-part exceeding E. Bino- 
mial confidence intervals are applicable to each point of the 
cumulative distribution as there are just two values of interest, 
e.g., satisfactory or unsatisfactory. P is a special case where 
E = 0. The robustness metric resulting from the cumulative 
probability distribution is directly related to classical concepts 
of rates of decay (growth) of first- and second-order closed- 
loop responses, time-to-half, and time-to-double. Taking de- 
gree-of-stability analysis further, rather than a vertical dis- 
criminant line, one can confine the closed-loop roots to sectors 
in the complex plane bounded by lines of constant damping 
and arcs of constant natural frequency. 6 Systems with roots 
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confined to these regions would be expected to display a cer- 
tain transient response speed. Again, the probability of roots 
lying within a sector follows a binomial distribution, and 
binomial confidence intervals apply. 

Performance specifications for aircraft flying qualities are 
detailed in Ref. 7 in terms of longitudinal and lateral-direc- 
tional criteria at three levels of performance for each flight 
phase. Many flying-qualities criteria require little computation 
above and beyond eigenvalue computation, making perfor- 
mance robustness as easy to characterize as stability robust- 
ness. For example, the short-period response can be character- 
ized by its damping ratio and natural frequency vs normal 
acceleration sensitivity to angle-of-attack n Q . The latter is 
illustrated in Ref. 7 by plotting the short-period undamped 
natural frequency vs n at as shown in Fig. 1. n a is simply a 
function of the dynamic pressure Q and vehicle parameters 


= 


^c L 

mg 


U) 


C Lm is the lift-curve slope, S ref the wing reference area, m the 
mass, and g the gravitational constant. Short-period-mode 
requirement levels for each flight phase are characterized by 
calculating the closed-loop eigenvalues and evaluating Eq. !. 
Repeated evaluations using Monte Carlo analysis give a distri- 
bution that can be shown pictorially on Fig. 1; the resulting 
measure of performance robustness is the probability of re- 



Flg.I Short-period response rs characterized by n a vs for cate- 
gory B flight phase (climb, cruise, descent) and mJi aircraft classes. 7 





Fig. 2 Example of step response bounds formed by scalar perfor- 
mance characteristics. 



Fig. 3 Confidence Interval calculation on the difference A P between 
two probabilities P\ and Pi. 


Table 1 Longitudinal parameters of tbe twin-jet aircraft 


Uniform 

variation* 

Description 

15 

Mass, slugs 

15 

Moment of inertia about they axis, slug- ft 2 

2 

Wing reference area, ft 

2 

Aerodynamic chord, ft 

2 

Wing span, ft 

30 

Center -of-gravity location as a percent of mean aerody- 
namic chord 

25 

Lift-curve slope 

25 

Lift-curve intercept 

40 

Deviation of the basic lift coefficient due to Mach effects 
on lift-curve intercept 

40 

Deviation of the basic lift coefficient due to Mach effects 
on lift-curve slope 

5 

Variation in lift coefficient with rate oT change of nondi- 
mensiona! a 

7.5 

Variation in lift coefficient with rate of change of nondi- 
mensiona! q 

10 

Variation in lift coefficient with change in elevator angle 

50 

Basic low-speed drag coefficient 

25 

Moment-curve slope 

25 

Moment-curve intercept 

25 

Deviation of the basic moment coefficient due to Mach 
effects on moment-curve intercept 

10 

Deviation in the basic moment coefficient due to Mach 
effects on moment-curve slope 

10 

Variation in moment coefficient with rate of change of 
nondimensional or 

10 

Variation in lift coefficient with rate of change of nondi- 
mensional q 

15 

Variation in moment coefficient with change in elevator 
angle 

10 

Cemer-of-gravity variation factor 


a ± percent of nominal parameter value 


maining within level 1, 2, or 3 criteria. 7 Binomial confidence 
interval computations can be applied to the scalar probability 
estimate. 

Time responses provide the most clear-cut means of evaluat- 
ing performance. Stochastic performance robustness can be 
portrayed as a distribution of possible trajectories around a 
nominal or desired trajectory. After defining “envelopes” 
around the nominal trajectory (Fig. 2), the probability of 
violating the envelopes can be computed using Monte Carlo 
evaluation. The envelope chosen around the nominal trajec- 
tory encompasses scalar performance measures; the trajecto- 
ries in Fig. 2 are examples of bounds defined by minimum 
and/or maximum allowable dead time, delay time, rise time, 
time-to-peak overshoot, peak overshoot, settling time, and 
steady-state error. 6 Although it is simple to conclude that a 
response violates an envelope, individual responses within the 
envelope may not be acceptable. In such cases, the derivative 
of a response and envelopes around the derivative also can be 
used as performance criteria. 3 

The criteria defining envelopes that bound an acceptable 
time response are not unique; the segmented envelopes in Fig. 
2 can be smoothed, or other scalars can be used to define 
points on the envelope. However, once an envelope is defined, 
time response distributions due to a command input, distur- 
bance, initial condition, or some combination can be com- 
puted by Monte Carlo methods. For each evaluation, the tra- 
je ctory is a binomial variable; it e ither st ays within the envel- 
ope or violates the envelope, and binomial con fidence inter- 
vals apply. Although individual time responses require more 
computation time than do individual sets of eigenvalues, such 
analysis is well within the capability of existing workstations. 

Confidence intervals for the difference between two proba- 
bilities are useful when comparing two control system designs. 
A statistic on the difference decides whether one controller is 
more robust than another, either as part of an iterative design 
process or as imbedded in an optimization technique. The 




Table 2 Scalar performance criteria defining 
command response envelope 


Scalar metric 

Value 

Maximum dead time 

2.5 s 

Maximum nonminimum- 

phase response 

-0.1 of desired steady-state value 

Minimum and maximum 

delay time 

1 .0 s and 7.5 s 

Minimum and maximum 

rise time 

2.0 s and 15.0 s 

Minimum and maximum 


peak time 

3.0 s and 18.0 s 

Maximum peak overshoot 

1 .25 of desired steady-state value 

Maximum settling time 

22.0 s 

Minimum and maximum 

steady-state error 

±0.025 of desired steady-state value 


Table 3 Setpoint for individual velocity 
and flight-path-anglr commands 


Command 

6T y *k 

deg 

y, fps 

7. deg 

q, rad/s 

a, deg 

V = 15 fps 

1.1 

15.3 

15 

0 

0 

-0.25 

7=4 deg 

24.1 

0.6 

0 

4 

0 

-0.01 


statistics literature gives several methods of computing the 
confidence interval for the difference between two binomial 
variables. Reference 8 presents a method based solely on indi- 
vidual confidence intervals. Given individual intervals based 
on independent Monte Carlo trials, 


Pr(I,<P, <£/,) = 1 — Qt| 

(2) 

Pr (l 2 <P 2 *U 2 )= 1 - or 2 

(3) 


the confidence interval around AP £ P x - P 2 is given by 8 

Pr[(T, - U 7 )<AP<{U l -L:)) 

2 I - aj - a 2 + (4) 

When identical parameter sets are used to generate individual 
intervals, the right-hand side of Eq. (4) is 1 - arj - cr 2 . Since 
(I,, U\) and (Z. 2 , C/ 2 ) are computed using the binomial test and 
represent exact intervals for the individual estimates, Eq. (4) is 
not an appioximation. Confidence interval comparisons are 
illustrated schematically in Fig. 3. The interpretation of the 
confidence interval for the difference is straightforward; the 
probability that the true difference lies within [(E, - l/ 2 ), (t/| 
- L 2 ) ] is at least 1 - o l - a 2 + a } a:. If the interval on A P 
contains zero {i.e., if the individual intervals overlap as they 
do initially in Fig. 3), then the difference in robustness be- 
tween the two systems is not proven significant at that number 
of evaluations. If the true difference AP is small, a larger 
number of evaluations may result in an interval that does not 
contain zero, as in Fig. 3. 

A given AP can result from many combinations of individ- 
ual probability estimates, and it is difficult to generalize the 
number of evaluations necessary to detect a difference of a 
certain magnitude. Nevertheless, the number of evaluations 
required for an individual confidence interval can be used to 
foretell the number of evaluations necessary to detect a differ- 
ence between two estimates. Figure 4 gives the required num- 
ber of evaluations J for each individual confidence interval, 
for the special case, a, = a 2 = 0.05. Using the difference 
P 2 - P, as the ordinate and P, as the abscissa, the curves show 
the minimum number of evaluations required to establish a 
significant difference. For example, if the probability esti- 
mates (denoted P) are P 2 - 0.45 and P\ - 0.4, Fig. 4 shows 
that a statistically significant difference (i.e., nonoverlapping 


confidence intervals) can be determined using approximately 
1500 Monte Carlo evaluations. Individual estimates of P 2 * 
0.15 and P x = 0.1 result in the same difference, but fewer than 
750 evaluations are required to detect the difference. Figure 4 
is based on individual confidence interval calculations, as 
presented in Ref. 3. 

Performance Robustness of Longitudinal Controllers 
For a Jet Transport 

SRA is applied to a twin-jet transport aircraft, with the goal 
of characterizing the performance robustness of longitudinal 



Fig. 4 Number of evaluations establishing significant differences be- 
tween two probabilities for 95^# confidence intervals and equal num- 
bers of evaluations for individual probabilities. 



a) Stochastic root locus with sector bounds defined by minimum 
level 1 short-period damping for cruise flight 



b) Short-period frequency vs acceleration sensitivity distribution 

Fig. 5 Stochastic robustness evaluation of the open-loop short-pe- 
riod dynamics of the twin-jet aircraft, based on 10,000 Monte Carlo 
evaluations. 


command responses. The rigid-body nonlinear longitudinal 
equations are 



- D + T cos(a) 
m 

L + T cos(of) 
mV 

M 

Iyy 

q - y 


- g sin(y) 

g cos(y) 
V 


(5) 


where [ V , y, q, a] represent velocity, flight-path-angle, pitch 
rate, and angle-of-attack, [L,D,M] are aerodynamic lift, 
drag, and pitching moment, T is the thrust, and g is the grav- 
itational constant. Equation (5) depends on a number of pa- 
rameters given in Table 1. Mean parameter values of the 
stability derivatives in Table l are functions of Mach number 
and altitude; they are interpolated from aerodynamic data 
curves for the aircraft at a given trim condition. 9 The aerody- 
namic model used to compute L, D, and M is a simplified 
version of that given in Ref. 9, modified to use only two lon- 
gitudinal controls (thrust and elevator). In this example, each 
Monte Carlo evaluation begins with the nonlinear equations 
of motion and associated parameters. The nonlinear equations 
are evaluated using appropriately distributed random parame- 
ters and are then linearized around the nominal trim condi- 
tion. The closed-loop eigenvalues and performance metrics are 
evaluated from the linearized system. 

The parameters are assumed to have uniform variations of 
the magnitudes given in Table I . For the wing parameters (5 ref , 
chord, span), these variations are representative of loose man- 
ufacturing tolerances. The mass and moment-of-inertia varia- 
tions are based on the maximum and minimum possible values 
of these parameters given in Ref. 9. The remaining parameter- 
variation estimates are based on interpolation accuracy and 
possible flight condition variations around the nominal value. 

25.0 


15.0 


> 

5.0 


time (sec) 

a) IS fps velocity commind: velocity response 



Trim conditions for a flight condition of V = 425 fps (130 
m/s) at an altitude of 5000 ft (1524 m) are as follows: thrust 
= 27.3^b, elevator = -0.65 deg, and angle-of-attack = 2.15 
deg. The open-loop eigenvalues for the state matrix resulting 
from linearizing Eq. (5) around trim are X = - 1.32 ± 2.44y, 
-0.0053 ± 0.0%2y. Stochastic robustness evaluation using 
the short-period Mil-spec requirements 7 shows an acceptable 
open-loop short-period mode for the uniform parameter vari- 
ations given in Table 1. Figure 5a shows the stochastic root 
locus with sectors defined by minimum level 1 short-period 
damping ratio for cruise or climb (category B flight phase); for 
10,000 evaluations, the short-period eigenvalues never violate 
the level 1 damping restriction. Figure 5b characterizes the 
short-period frequency vs acceleration sensitivity, which also 
remains within level 1 constraints for 10,000 evaluations. The 
probability estimate of violating level 1 short-period specifica- 
tions is 0, with 95^o confidence intervals of (0, 3.69E - 4). 

Design of Longitudinal Controllers 
A command response that stays within the envelope de- 
scribed T>y scalar criteria in Table 2 serves as the performance 
requirement for designing linear regulators for velocity and 
flight-path-angle commands. In addition, elevator deflections 
are limited to ± 30 deg, and thrust commands must remain 
between 0 and lOO^o. The desired commands y* - V* or 
y* = y* and corresponding setpoints x* - [V y q a] r , u* * 
(5 7” 5E] are given in Table 3. The open-loop responses to in- 
dividual velocity and flight-path-angle commands are inade- 
quate because of the slow, lightly damped phugoid mode. 
Numerical values ofjhe results that follow depend heavily on 
the performance criteria chosen. The envelopes defined In 
Table 2 reflect tolerable variations around an acceptable nom- 
inal response. The control limits are typical of those for a jet 
transport. Changing the time response envelopes or control 
authority limits would give different numerical results. The 
emphasis in this example is not on the specific criteria chosen, 
but on how SRA characterizes performance given a control 
system design and performance specifications. 



c) 15 fps velocity command: elevator response 



b) 15 fps velocity command: thrust response d) 4-deg flight-path-angle command: flight-path-angle response 

Fig. 6 Closed-loop command responses using IMF controller, 500 Monte Carlo evaluations. Nominal response Is indicated by the solid line. 
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Structured linear-quadratic regulators 10 offer a simple means 
of designing a linear control system with desirable perfor- 
mance and robustness characteristics. Specifications of the lin- 
ear-quadratic performance index and subsequent control gains 
using implicii-model-following (IMF) minimizes the dynamic 
response error between the closed-loop system and an ideal 
model. 10 State, control, and cross-weighting matrices (Q, R, 
M) are^ased^on a quadratic cost function that weights the 
difference between the actual state rate (x) and that of an ideal 
model (x M ), where 

x M = F w x M (6) 


IMF offers a straightforward way of designing controllers that 
approximate desired dynamic characteristics. For this exam- 
ple, the ideal model was chosen to increase the natural fre- 
quency and damping of the phugoid mode, while maintaining 


acceptable short period response: 


-0.3 

-32.17 

-0.0104 

-23.34 

0.00381 

-0.1949 

0.0006 

1.356 

0.0 

-0.0 

-1.273 

-5.981 

0.0038 

0.1949 

0.999 

-1.356 


f - 1-35 

± 2.39 j 


7m = 

(-0.213 

±0.314y 



(7) 


( 8 ) 



time (sec) 


a) 15 fps velocity command: velocity response 



b) 15 fps velocil) command: thrust response 


Stochastic performance robustness analysis is based on the 
probability of violating the desired time response envelopes 
(A and A) and the probability of control saturation (Ar and 

The IMF controller gives a nominal closed-loop command 
response to separate velocity (Figs. 6a-c) and Right -path-angle 
(Fig. 6d) commands that is within the acceptable time-re- 
sponse envelope. Figure 6 also shows 500 Monte Carlo evalua- 
tions of the command response; the nominal steady-state con- 
trol inputs and state arc given in Table 3, and the nominal 
response in Fig. 6 is indicated by a solid line. The response and 
associated envelopes in Fig. 6 are shown for the commanded 
variable only; the remaining state elements do not require 
performance constraints in this example. Thrust and elevator 
time histories are shown for the velocity command response 
only. Parameter uncertainty effects appear as variations 
around the nominal response, indicated by the dark distribu- 
tion and associated outliers. Parameter uncertainty results in a 
distribution of transient responses that stays within the envel- 
ope, and nonzero steady-state errors that violate the envelope 
for both velocity (Fig. 6a) and flight-path-angle (Fig. 6d) 
commands. Based on 500 Monte Carlo time response evalua- 
tions, the estimate A is 0.002 with 95% confidence intervals 
(5. IE - 5,0.0111) and the estimate A is 0.368 (0.326,0.412). 
The nominal elevator response violates control limits for both 
command responses, and in each case, the probability of 
elevator saturation is Az - 1.0. Note that the control satura- 
tion limits in Figs. 6b-c are adjusted to reflect the remaining 
control authority after considering trim requirements. 



time (sec) 


d) 4-deg fllght-path-angle command: fllght-path-angle response 



e) 4-deg flighl-path-angle command: thrust response 



time (sec) time (sec) 

c) 15 fps velocity command: elevator response f) 4-deg flighl-patb-angle command: elevator response 


Fig. 7 Closed-loop command response using PF1MF controller, with filler control weighting R f * diagOO, 50), 500 Monte Carlo evaluations. 
Nominal response Is Indicated by the solid line. 
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Fig . S Stochastic performance robustness evaluation with PFIMF: 
Probability of violating flight-path-angle command response A T and 
probability of violating elevator saturation limits Pie vs filler weight 
Solid lines give probability estimates, dashed lines give confi- 
dence intervals. 
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time (sec) 

Fl|.* Closed-loop command response using PFIMF controller, 
with niter control weighting Rjr =diag(10, 50), 500 Monte Carlo eval- 
uations: 15 fps velocity command subject to constant disturbance 
w 9 * 40 fps. Nominal response is indicated by the solid line. 



Implicit model following modified by state augmentation 10 
can help meet control authority constraints. Proportional-fil- 
ter (PF) compensation adds integrators to restrict the control 
rates, thus preventing instantaneous control changes and re- 
ducing the maximum control effort. The control vector is 
appended to the state vector 



where F and G are the nominal dynamic and control effect 
matrices, % = x{t)~x*, A = u(f) - u # , and v(0 is a com- 
manded control rate. The PFIMF state weighting matrix is 


Qf = 


Q 

M r 


M 

R 


00 ) 


where Q, R, M are the original (IMF) weighting matrices. A 
weighting matrix, R £ , constrains the control rates. Elements 
of R/- affect the bandwidth of each control; the larger the 
weight, the more the control rate is restricted. 

The IMF regulator is augmented to include low-pass filter- 
ing of the control command, with a diagonal control-rate 
weighting matrix R F =* diag[10, 50). Figure 7 shows 500 
stochastic state and control histories to individual velocity and 
flight-path-angle commands using the PFIMF controller and a 
stream of random numbers independent from the IMF case. 
The (1, 1) element of R f (fl 4r ) determines the amount of 
filtering on thrust rate, and the (2, 2) element (/? 4£ ) controls 
elevator rate. With filter elements, the control rates are no 
longer unlimited, and the mean control responses remain un- 
saturated. Steady-state error due to parameter uncertainty 
remains within the desired state history envelope for the veloc- 
ity command response (Fig. 7a). Steady-state error for the y 


command improves, although the variation in the y transient 
response is much greater than that of the IMF regulator alone, 
as seen by comparing Figs. 6d and 7d. P v and P y estimates 
corresponding to Fig. 7 are 0.0 (0.0, 0.0074) and 0.034 (0.0199, 
0.0539), respectively. For 500 evaluations, the PFIMF flight- 
path-angle command response improvement over the IMF 
case alone proves significant by application of confidence in- 
tervals on the difference (P^imf - f’yptMF)- Applying Eq. 4, 


Pr [0.2721 s (Py, MF - Ptpfimf) * 0.3921] a 0.9025 (1 1) 

Equation 1 1 states that with PF augmentation between 27 and 
39^«, more of the flight-path-angle responses lie within the 
envelope, with a confidence coefficient of at least 0.9025. The 
mean elevator response for the flight-path-angle command 
dips just to saturation limits, and the probability of elevator 
saturation is = 0.502 (0.457, 0.547). 

Stochastic robustness analysis shows that PF augmentation 
improves performance objectives by reducing control rates 
and steady-state error due to uncertainty. The state and con- 
trol response to the velocity command prove acceptable {P v> 
P 4£ , and \ Pij all equal 0), and the improved responses to flight- 
path-angle command are statistically significant. For the 
flight-path-angle command, SRA demonstrates the tradeoff 
between the two performance objectives; increasing the (2, 2) 
element (/? 4£ ) of R £ will further reduce elevator command 
authority at the expense of the y time response. Figure 8 
illustrates this tradeoff by showing P y , and their confl- 
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Fig. 10 Closed-loop command response using PIF1MF controller, 
with filter control weighting Rp * dlog(200, 50), and integral stale 
weighting Q/ = diag(0.1, 100), 500 Monte Carlo evaluations: 15 fps 
velocity command subject to constant disturbance w ¥ ■ 40 fps. Nomi- 
nal response is indicated by the solid line. 
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Fit* II Stochastic performance robustness evaluation with PIF1MF: 
Probability of violating velocity command response envelope Pv and 
probability of violating thrust saturation limits P iT vs filter weight 
Rij, Solid lines give probability estimates, dashed lines give confi- 
dence Intervals. 


dcnce intervals as functions of the design parameter R A 
plot like Fig. 8 can be used to choose the filter weight that 
gives the smallest probabilities of envelope violation while 
adhering as well as possible to the control authority restric- 
tions. In this case, it is not possible to simultaneously reduce 
P y and Pie to zero by varying R 6£ . Nevertheless, stochastic 
robustness analysis offers a simple, understandable means of 
relating design parameters to performance objectives and of 
choosing the best control gains to meet those objectives. 

Design of a Longitudinal Controller 
for Disturbance Rejection 

As a final example, the preceding analysis is extended to 
encompass a performance constraint on disturbance rejection. 
The equations of motion are modified to include a vertical 
wind disturbance w r 

r-D sin{cr - a„) + T cos(or - a*) . , " 

i g sm(y) 

m (12) 

L sin (or - a g ) + T COS(a - ct a ) _ g CQS(y) 

mV V 



where 


a 0 


a 4 y — tan “ 1 


V sin(>) + u\ 
V cos(-y) 


(13) 


With the disturbance present, the state components represent 
inertial velocity, flight-path-angle, pitch, and angle-of-attack, 
and the disturbance enters through the expression for air-rela- 
tive angle-of-attack a fl . A disturbance input matrix is defined 
for robustness analysis by numerical linearization of the non- 
linear equations with respect to around the nominal condi- 
tion w v =0. Velocity command response subject to a constant 
40-fps vertical velocity disturbance using the PFIMF controller 
is shown in Fig. 9. The mean response shows a nonzero steady- 
state error that violates the command response envelope, and 
uncertainty causes a larger spread around the nominal re- 
sponse than that of the system without the disturbance (Fig. 7). 
Also, the steady-state flight-path-angle (not shown) is less than 
zero due to the disturbance. 

Proportional-integral (PI) compensation introduces a com- 
mand-error integral for each commanded state element, zero- 
ing steady-state error and improving disturbance rejection 
characteristics. The perturbation equations for the nominal 
system are 


F 0 
H 0 


H 



0 0 o" 

1 0 0 . 


(14) 

(15) 


*</) = {(0) + j f(r)dr (16) 

and y(f) = y(/)- y*. Here, y* * \V y] r , and a (2 x 2) weighting 
matrix Q/ is appended to the original state weighting matrix. 
Diagonal elements of Q/ affect the rate at which the command 
error integrals approach zero. The diagonal components are 
chosen to keep the velocity command within the desired envel- 
ope and to zero the flight-path-angle response. Command er- 
ror integrals are added to the existing PFIMF controller, and 
for the resulting PIFIMF system with Q, - diag[0.01 , 1001 and 
R f = diag[200, 50], Fig. 10 shows an improved velocity com- 
mand response y* = (P*0] r . The 500-cvaluation probability 
estimates and 95V« confidence intervals are Py *0 (0.0, 7.4E- 
3) and P hT = 0.002 (5.1E-5, 0.0111). The (1, I) component of 
R F is increased to restrain thrust as the command error inte- 
grals are introduced. Figure 1 1 shows analysis of the tradeoff 
between Py and P hT a function of design parameter R*r 
comparable to that presented for the flight-path-angle re- 
sponse in Fig. 8. Again, Fig. 1 1 can be used to choose control 
system design parameters that best meet performance objec- 
tives. 


Conclusion 

Stochastic robustness analysis offers a rigorous yet straight- 
forward alternative to other robustness metrics that is simple 
to compute and is unfettered by normally difficult problem 
statements, such as non-Gaussian statistics, products of pa- 
rameter variations, and structured uncertainty. The analysis 
embraces both stability and performance metrics, handling 
qualities requirements, and more general responses. Binomial 
confidence intervals provide statistical bounds on the proba- 
bility of instability and on performance metrics. Statistical 
comparisons of control system robustness also are rendered 
through confidence intervals. Both stability and performance 
metrics resulting from stochastic robustness analysis provide 
details relating system specifications intrinsic to a given appli- 
cation and control system design parameters. Stochastic ro- 
bustness analysis has a significant role to play in computer- 
aided control system design. 
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